GÉO-VISUALISATION AVEC R
Action nationale de formation
MATE SHS
GÉO-VISUALISATION AVEC R
Action nationale de formation
MATE SHS
Cette formation de 3 heures 30 porte sur la visualisation de données géographiques sous R. Y sont abordées les traitements SIG de base, la cartographie thématique (cartes en figurés proportionnels, cartes choroplèthes, cartes de typologie, etc) et des cartes reposant sur des techniques plus avancées comme les cartes sur grille ou les cartes de discontinuité.
Sont ici principalement utilisés les packages sf (manipulation de données spatiales) et cartography (cartographie thématique). Chaque exemple propose une chaine de traitement depuis le chargement des données, leur mise en forme, jusqu’à la mobilisation de méthodes adaptées permettant de répondre à des questions spatiales.
Les exemples proposés sont basés sur différentes échelles, du local (région Occitanie) en passant par l’Europe et le Monde.
Packages & données
Environnement de travail
Pour exécuter les programmes proposés par la formation, installez / mettez à jour R (version 3.4.4 minimum), R Studio. Lancez ensuite R studio et installez les packages au moyen des commandes suivantes (prévoir bien 10 minutes)
Les packages
Installation des packages utilisés pendant la formation.
install.packages("sf")
install.packages("cartography")
install.packages("countrycode")
install.packages("readxl")
install.packages("rnaturalearth")
install.packages("cartogram")
install.packages("reshape2")
install.packages("eurostat")Optionnel
install.packages("rmapshaper")
install.packages("osmdata")Version des packages utilisés dans ce document
## [1] "sf (version 0.6.3), cartography (version 2.1.2), countrycode (version 1.0.0), readxl (version 1.1.0), rnaturalearth (version 0.1.0), cartogram (version 0.1.0), reshape2 (version 1.4.3), eurostat (version 3.2.9), rmapshaper (version 0.4.0), osmdata (version 0.0.7.1)"
Packages dédiés à visualisation de données géographiques
- Le package cartography permet de créer et intégrer des cartes thématiques dans sa chaine de traitement en R. Il permet des représentations cartographiques tels que les cartes en symboles proportionnels, des cartes choroplèthes, des typologies, des cartes de flux ou des cartes de discontinuité. Il offre également des fonctions qui permettent d’améliorer la réalisation de la carte, comme des palettes de couleur, des éléments d’habillage (échelle, flèche du nord, titre, légende…), d’y rattacher des labels ou d’accéder à des APIs cartographiques.
- Le package cartogram comme son nom l’indiqie permet de construire des cartogrames continus ou non continus.
Packages dédiés à la manipulation de données (spatiales ou non)
- Le package sf est un package qui permet d’importer, gérer et transformer des données géographiques (gestion des projections, opérations SIG).
- Le package readxl permet d’importer des fichiers Excel sous R.
- Le package reshape2 permet de transformer et mettre en forme des données sous R.
- Le package rmapshaper permet de généraliser le contours d’un fond de carte
Packages fournisseurs de données
- Le package eurostat permet d’importer des données d’Eurostat, fournisseur de référence pour les territoires européens.
- Le package rnaturalearth permet d’importer les fonds de carte vectoriels de référence pour le Monde (niveaux État ou infra-nationaux) ainsi que des éléments d’habillage (graticules, traits de côte, etc.)
- Le package osmdata permet d’importer des données d’OpenStreetMap.
Téléchargez les données
Téléchargez les données sur la page github. Décompressez les données sur votre espace de travail.
Exo 1 - Bienvenue à Sète
Commandes de base
Définissez votre répertoire de travail.
setwd("Votre_repertoire_de_travail")Consulter le contenu. Et regarder ce qu’il y a dans le répertoire “data”
list.files()## [1] "data" "data.zip" "index.html" "index.Rmd" "outputs"
## [6] "pics" "prg"
list.files("data")## [1] "download" "Europe" "France" "Hérault" "Occitanie" "world"
Import d’une couche SIG dans R avec le package sf. Pacakge basé sur GEOS (Geometry Engine Open Source) et GDAL (Geospatial Data Abstraction Library).
library("sf")## Linking to GEOS 3.6.2, GDAL 2.2.3, proj.4 4.9.3
communes <- st_read(dsn = "data/Hérault/communes.shp", stringsAsFactors = F)## Reading layer `communes' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/Hérault/communes.shp' using driver `ESRI Shapefile'
## Simple feature collection with 343 features and 5 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 662675 ymin: 6234874 xmax: 796333 ymax: 6319348
## epsg (SRID): NA
## proj4string: +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
Voir la table atributaire
communes## Simple feature collection with 343 features and 5 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 662675 ymin: 6234874 xmax: 796333 ymax: 6319348
## epsg (SRID): NA
## proj4string: +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
## First 10 features:
## INSEE_COM NOM_COMM STATUT SUPERFICIE
## 1 34029 BELARGA Commune simple 415
## 2 34290 SAINT-VINCENT-DE-BARBEYRARGUES Commune simple 222
## 3 34032 BEZIERS Sous-prfecture 9567
## 4 34082 COMBAILLAUX Commune simple 902
## 5 34108 FRONTIGNAN Chef-lieu canton 3984
## 6 34166 MONTBLANC Commune simple 2716
## 7 34066 CAZEVIEILLE Commune simple 1619
## 8 34296 SAUSSINES Commune simple 630
## 9 34269 SAINT-JEAN-DE-MINERVOIS Commune simple 3277
## 10 34012 ARGELLIERS Commune simple 5065
## NOM_DEPT geometry
## 1 HERAULT POLYGON ((742075 6272005, 7...
## 2 HERAULT POLYGON ((770896 6289405, 7...
## 3 HERAULT POLYGON ((718759 6244261, 7...
## 4 HERAULT POLYGON ((761606 6284489, 7...
## 5 HERAULT POLYGON ((758738 6257679, 7...
## 6 HERAULT POLYGON ((728067 6248191, 7...
## 7 HERAULT POLYGON ((762395 6294296, 7...
## 8 HERAULT POLYGON ((784582 6294234, 7...
## 9 HERAULT POLYGON ((686897 6252592, 6...
## 10 HERAULT POLYGON ((756922 6287661, 7...
head(communes)## Simple feature collection with 6 features and 5 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 710431 ymin: 6244261 xmax: 771300 ymax: 6292043
## epsg (SRID): NA
## proj4string: +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
## INSEE_COM NOM_COMM STATUT SUPERFICIE
## 1 34029 BELARGA Commune simple 415
## 2 34290 SAINT-VINCENT-DE-BARBEYRARGUES Commune simple 222
## 3 34032 BEZIERS Sous-prfecture 9567
## 4 34082 COMBAILLAUX Commune simple 902
## 5 34108 FRONTIGNAN Chef-lieu canton 3984
## 6 34166 MONTBLANC Commune simple 2716
## NOM_DEPT geometry
## 1 HERAULT POLYGON ((742075 6272005, 7...
## 2 HERAULT POLYGON ((770896 6289405, 7...
## 3 HERAULT POLYGON ((718759 6244261, 7...
## 4 HERAULT POLYGON ((761606 6284489, 7...
## 5 HERAULT POLYGON ((758738 6257679, 7...
## 6 HERAULT POLYGON ((728067 6248191, 7...
View(communes)L’instruction plot permet d’afficher la couche. L’instruction st_geometry permet d’accéder à la variable définissant les géométries.
plot(st_geometry(communes))Afficher la couche avec des parametres graphiques
plot(st_geometry(communes), col="#aec8f2", border="darkblue", lwd=1)Introduction au SIG avec R / Manipuler les informations spatiales
Nous cherchons à identifier les communes qui ont une partie de leur territoire situé à moins de 20 km du centre de Sete.
Nous commençons par extraire la commune de Sète.
macommune <- "SETE"
monpoly <- communes[communes$NOM_COMM == macommune,]
plot(st_geometry(communes), col="#aec8f2", border="darkblue", lwd=1)
plot(st_geometry(monpoly), col="red", border="purple", lwd=1, add=T)Extraitre le centroide de la commune de Sete.
moncentre <- st_centroid(x = monpoly)## Warning in st_centroid.sf(x = monpoly): st_centroid assumes attributes are
## constant over geometries of x
plot(st_geometry(communes), col="#aec8f2", border="darkblue", lwd=1)
plot(st_geometry(monpoly), col="red", border="purple", lwd=1, add=T)
plot(st_geometry(moncentre), pch=20, col="black", cex=2, add=T)Calculer une zone tampon
mydist <- 20000
buff <- st_buffer(x = st_geometry(moncentre), dist=mydist)
plot(st_geometry(communes), col="#aec8f2", border="darkblue", lwd=1)
plot(st_geometry(monpoly), col="red", border="purple", lwd=1, add=T)
plot(st_geometry(moncentre), pch=20, col="black", cex=2, add=T)
plot(buff, col=NA, border="black", lty=2,add=T)Récupérer la liste des communes
communes$buff <- st_intersects(st_geometry(communes), st_geometry(buff), sparse = FALSE)
head(communes)## Simple feature collection with 6 features and 6 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 710431 ymin: 6244261 xmax: 771300 ymax: 6292043
## epsg (SRID): NA
## proj4string: +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
## INSEE_COM NOM_COMM STATUT SUPERFICIE
## 1 34029 BELARGA Commune simple 415
## 2 34290 SAINT-VINCENT-DE-BARBEYRARGUES Commune simple 222
## 3 34032 BEZIERS Sous-prfecture 9567
## 4 34082 COMBAILLAUX Commune simple 902
## 5 34108 FRONTIGNAN Chef-lieu canton 3984
## 6 34166 MONTBLANC Commune simple 2716
## NOM_DEPT geometry buff
## 1 HERAULT POLYGON ((742075 6272005, 7... TRUE
## 2 HERAULT POLYGON ((770896 6289405, 7... FALSE
## 3 HERAULT POLYGON ((718759 6244261, 7... FALSE
## 4 HERAULT POLYGON ((761606 6284489, 7... FALSE
## 5 HERAULT POLYGON ((758738 6257679, 7... TRUE
## 6 HERAULT POLYGON ((728067 6248191, 7... FALSE
Récupérer le nombre de communes à moins de 20 km
nb <- dim(communes[communes$buff == T,])
nb[1]## [1] 40
Extraction et affichage des communes
communes20km <- communes[communes$buff == TRUE,]
plot(st_geometry(communes), col="#aec8f2", border="darkblue", lwd=1)
plot(st_geometry(monpoly), col="red", border="purple", lwd=1, add=T)
plot(st_geometry(moncentre), pch=20, col="black", cex=2, add=T)
plot(buff, col=NA, border="black", lty=2,add=T)
plot(st_geometry(communes20km), col=NA, lwd=2, border="red", add=T)Ajout d’un titre
plot(st_geometry(communes), col="#aec8f2", border="darkblue", lwd=1)
plot(st_geometry(monpoly), col="red", border="purple", lwd=1, add=T)
plot(st_geometry(moncentre), pch=20, col="black", cex=2, add=T)
plot(buff, col=NA, border="black", lty=2,add=T)
plot(st_geometry(communes20km), col=NA, lwd=2, border="red", add=T)
montitre <- paste0("Il y a ", nb[1], " communes situées à moins de ", mydist/1000, "km de ",tolower(macommune))
title(montitre)Export (si besoin)
st_write(obj = communes20km, dsn = "outputs/resultat.shp")A vous de jouer !
Reproduisez cette procédure avec la commune et la distance de votre choix.
Pour accéder au listing des communes d’Hérault
communes$NOM_COMMExo 2 - Cartographier l’Occitanie
Objectifs
- Thématique : Réaliser une carte des catégories Socio professionnelles.
- Espace d’étude : Occitanie, communes
- Visualisation associée : carte en figurés proportionnels (+ cartogrammes de Dorling)
- packages utilisés : sf, cartography, cartogram, readxl
Créer un nouveau fichier “Exo2.R”.
Chargement et préparation des données
Pour bien commencer, il est fondamental de jeter un oeil au fichier d’entrée pour comprendre son organisation et savoir comment l’importer sous R.
Ouvrez le fichier data/France/INSEE/pop-act2554-csp-cd-6814.xls Nous nous intéressons ici à l’onglet “COM_2014”.
library("sf")
library("readxl")Importez correctement cette feuille du tableau Excel dans R. Les commandes de base si dessous permettent de visualiser le nom des colonnes, de filtrer le contenu de la table (par région) et de visualiser finalement le tableau de données dans R.
sheet <- "COM_2014"
data <- data.frame(read_excel("data/France/INSEE/pop-act2554-csp-cd-6814.xlsx", skip = 15, sheet = sheet))
colnames(data)## [1] "RR" "DR"
## [3] "CR" "STABLE"
## [5] "DR16" "LIBELLE"
## [7] "csx_rec1taxtypac_rec1rpop2014" "csx_rec1taxtypac_rec2rpop2014"
## [9] "csx_rec2taxtypac_rec1rpop2014" "csx_rec2taxtypac_rec2rpop2014"
## [11] "csx_rec3taxtypac_rec1rpop2014" "csx_rec3taxtypac_rec2rpop2014"
## [13] "csx_rec4taxtypac_rec1rpop2014" "csx_rec4taxtypac_rec2rpop2014"
## [15] "csx_rec5taxtypac_rec1rpop2014" "csx_rec5taxtypac_rec2rpop2014"
## [17] "csx_rec6taxtypac_rec1rpop2014" "csx_rec6taxtypac_rec2rpop2014"
data <- data[data$RR == "76",] # Ne prendre que les communes d'Occitanie (code 76)
View(data) # vOIR Si le fichier a correctement été importé. Pour notre analyse nous avons besoin de regrouper plusieurs colonnes. Créez de nouveaux champs en aggrégeant emploi + chômage par CSP. Renommez ces champs de façon explicite. Cela facilitera leur traitement par la suite. Agrégez également le code départemental au code communal afin de créer un identifiant (id) compatible avec le fond de carte IGN.
data$id <- paste0(data$DR, data$CR) # Creation du code commune
data[,"agr"]<- round(data[7] + data[8],0) # Actifs + chomeurs par CSP
data[,"art"] <- round(data[9] + data[10],0)
data[,"sup"] <- round(data[11] + data[12],0)
data[,"int"] <- round(data[13] + data[14],0)
data[,"emp"] <- round(data[15] + data[16],0)
data[,"ouv"] <- round(data[17] + data[18],0)
data[,"ouv_chom"] <- round(data[18],0)
data$total <- data$agr + data$art + data$sup + data$int + data$emp + data$ouvImportez et visualisez le fond de carte (extrait IGN Geofla communes).
communes <- st_read(dsn = "data/Occitanie/occitanie.shp", stringsAsFactors = F)## Reading layer `occitanie' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/Occitanie/occitanie.shp' using driver `ESRI Shapefile'
## Simple feature collection with 4516 features and 17 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 428145.7 ymin: 6137116 xmax: 848019.7 ymax: 6439628
## epsg (SRID): NA
## proj4string: +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
plot(st_geometry(communes),col="white", border="black", lwd=0.1)La Jointure est une étape cruciale qui permet de lier un fond de carte avec des données par un même identifiant. Il faut que les identifiants entre ces deux ressources soient exactement les mêmes.
Réalisez la jointure entre données et fond de carte avec l’instruction merge. Cela présuppose en amont que les identifiants géographiques entre ces deux ressources sont les mêmes (ça tombe bien, c’est le cas !). Identifiez bien le nom des champs sur lesquel doit porter la jointure.
communes <- merge(x = communes, y = data, by.x = "INSEE_COM", by.y = "id", all.x=T)Pour finir, nettoyez le tableau. On ne garde que les colonnes utiles pour la suite…
colnames(communes)## [1] "INSEE_COM" "ID_GEOFLA"
## [3] "CODE_COM" "NOM_COM"
## [5] "STATUT" "X_CHF_LIEU"
## [7] "Y_CHF_LIEU" "X_CENTROID"
## [9] "Y_CENTROID" "Z_MOYEN"
## [11] "SUPERFICIE" "POPULATION"
## [13] "CODE_ARR" "CODE_DEPT"
## [15] "NOM_DEPT" "CODE_REG"
## [17] "NOM_REG" "RR"
## [19] "DR" "CR"
## [21] "STABLE" "DR16"
## [23] "LIBELLE" "csx_rec1taxtypac_rec1rpop2014"
## [25] "csx_rec1taxtypac_rec2rpop2014" "csx_rec2taxtypac_rec1rpop2014"
## [27] "csx_rec2taxtypac_rec2rpop2014" "csx_rec3taxtypac_rec1rpop2014"
## [29] "csx_rec3taxtypac_rec2rpop2014" "csx_rec4taxtypac_rec1rpop2014"
## [31] "csx_rec4taxtypac_rec2rpop2014" "csx_rec5taxtypac_rec1rpop2014"
## [33] "csx_rec5taxtypac_rec2rpop2014" "csx_rec6taxtypac_rec1rpop2014"
## [35] "csx_rec6taxtypac_rec2rpop2014" "agr"
## [37] "art" "sup"
## [39] "int" "emp"
## [41] "ouv" "ouv_chom"
## [43] "total" "geometry"
fields <- c("INSEE_COM", "LIBELLE", "agr", "art", "sup", "int", "emp", "ouv","ouv_chom", "total", "geometry")
communes <- communes[,fields]
colnames(communes) <- c("id", "name", "agr", "art", "sup", "int", "emp", "ouv","ouv_chom", "total", "geometry")
head(communes)## Simple feature collection with 6 features and 10 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 556068.9 ymin: 6180702 xmax: 612006.7 ymax: 6219884
## epsg (SRID): NA
## proj4string: +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
## id name agr art sup int emp ouv ouv_chom total
## 1 09001 Aigues-Juntes 5 5 5 9 5 5 0 34
## 2 09002 Aigues-Vives 0 19 24 58 73 63 15 237
## 3 09003 Aiguillon 0 0 5 19 44 44 10 112
## 4 09004 Albiès 5 11 11 11 11 16 11 65
## 5 09005 Aleu 4 0 0 4 8 4 4 20
## 6 09006 Alliat 0 0 0 10 5 5 0 20
## geometry
## 1 MULTIPOLYGON (((572559.1 62...
## 2 MULTIPOLYGON (((607979.2 62...
## 3 MULTIPOLYGON (((612006.7 62...
## 4 MULTIPOLYGON (((594651 6188...
## 5 MULTIPOLYGON (((557843.5 61...
## 6 MULTIPOLYGON (((584509.3 61...
Visualisation : carte en figuré proportionnel
Le package cartography
Chargez le package cartography.
library("cartography")Ce package permet de créer et intégrer des cartes thématiques dans sa chaine de traitement en R. Il permet des représentations cartographiques tels que les cartes en symboles proportionnels, des cartes choroplèthes, des typologies, des cartes de flux ou des cartes de discontinuité. Il offre également des fonctions qui permettent d’améliorer la réalisation de la carte, comme des palettes de couleur, des éléments d’habillage (échelle, flèche du nord, titre, légende…), d’y rattacher des labels ou d’accéder à des APIs cartographiques.
Pour utiliser aisément ce package, plusieurs ressources d’intérêts
- La documentation du package qui documente toutes les fonctions du package, accessible directement dans R Studio. Pour celà, vous pouvez taper simplement :
?cartography- La cheat sheet de cartography, qui résume les principales fonctions du package de façon synthétique.
- La vignette associée au package, qui présente des réalisations issues de ce package, elle aussi accessible directement dans R studio.
- La vignette Faire des cartes - Introduction au package sf qui présente intelligemment les packages sf et cartography à partir d’un cas d’utilisation centré sur Mexico.
- Le blog R Géomatique, maintenu par l’auteur de cartography qui met à disposition ressources et exemples d’intérêt liés au package et à la représentation d’information spatiale sous R.
Éléments théoriques
Savoir représenter correctement de l’information dans des cartes thématiques, c’est d’abord savoir caractériser la nature statistique des données.
Jacques Bertin a notamment théorisé dans son ouvrage de référence “Sémiogie Graphique” (1967) les méthodes visuelles qui permettent de retranscrire efficacement une information statistique sous forme cartographique.
Il identifie des variables visuelles (ou rétiniennes) dont l’utilisation dépend de l’implantation (point, ligne, surface) et la nature statistique de l’information à représenter (association, sélection, ordre, quantité).
Pour les caractères quantitatifs absolus la variable visuelle adaptée est la taille afin de préserver graphiquement les rapports de proportionnalité existant entre les différentes valeurs du tableau de données.
Visualisation !
Cartographiez le nombre d’actifs par commune grâce à la fonction propSymbolsLayer
plot(st_geometry(communes), col="lightblue",border="white", lwd=0.2) # Fond de carte de la région Occitanie
propSymbolsLayer(x = communes, var = "total", # Carte en figuré proportionnel
symbols = "circle", col = "red",
legend.pos = "right", border = "grey",
legend.title.txt = "Les actifs",
legend.style = "c")
top <- sort(data.frame(communes)[,"total"], decreasing = TRUE)
labelLayer(x = communes[data.frame(communes)[,"total"] %in% top[1:15],], txt = "name", halo=TRUE, cex = 0.6, col= "#000000", bg = "#FFFFFF50", overlap = FALSE) # Nommer les 15 premières villes d'Occitanie en nombre d'actifs
title("Région Occitanie")Visualisation : carte combinant figuré proportionnel et typologie
Il est possible d’associer à ces figurés proportionnels la représentation d’un caractère qualitatif de différenciation.
La variable visuelle adaptée pour ce type de données est la couleur différencielle.
Exécutez le code ci-dessous, qui propose le calcul d’une typologie sur la surreprésentation d’ouvriers ou de cadres. Cette typologie peut être représentée simplement avec la fonction typoLayer.
Pour le choix des couleurs, vous pouvez utiliser color picker pour obtenir les couleurs que vous désirez au format hexadécimal présenté ci-dessous.
# Typologie (plus d'ouvriers que de cadres ?)
communes$typo <- "Indéterminé"
communes$r <- communes$ouv / communes$sup
communes$r[is.na(communes$r)] <- 0
communes[communes$r > 1.1,"typo"] <- "Plus d'ouvriers que de cadres"
communes[communes$r < 0.91,"typo"] <- "Plus de cadres que d'ouvriers"
# Choix des couleurs
colouv <- "#dd4e44" # tonalité rouge
colsup <- "#63b269" # tonalité verte
# Cartographie
typoLayer(communes, var="typo",
legend.values.order = c("Plus d'ouvriers que de cadres",
"Plus de cadres que d'ouvriers",
"Indéterminé" ),
col = c(colouv, colsup, 'grey'),
border = NA, legend.title.txt = "Type dominant") On peut aussi décider d’associer cette représentation à la carte précédente en figuré proportionnel grâce à la fonction propSymbolsTypoLayer :
# Combiner avec la carte précédente
plot(st_geometry(communes), col="lightblue",border=NA)
propSymbolsTypoLayer(x = communes, var = "total", var2 = "typo",
symbols = "circle",
inches = 0.4,
col = c(colouv, colsup, 'grey'),
border = "white", lwd = 0.1,
legend.var2.values.order = c("Plus d'ouvriers que de cadres",
"Plus de cadres que d'ouvriers",
"Indéterminé" ),
legend.var.pos = "right",
legend.var.title.txt = "Nombre d'actifs",
legend.var2.title.txt = "Type dominant")
labelLayer(x = communes[data.frame(communes)[,"total"] %in% top[1:15],], txt = "name", halo=TRUE, cex = 0.6, col= "#000000", bg = "#FFFFFF50", overlap = FALSE)
title("Cadres et ouvriers en région Occitanie")Visualisation : cartogramme de Dorling
A cette échelle géographique on remarque que les cercles se surimposent les uns aux autres. Cela ne nuit pas à la lisibilité générale de la carte mais on peut aussi décider d’utiliser un package et une fonction qui permet d’éviter cela.
Pour ce faire, lancez la librairie cartogram et la fonction cartogram_dorling.
library("cartogram")# cartogram de Dorling
actifs <- cartogram_dorling(communes, "total", k = 1, m_weight = 1, itermax = 20)
plot(st_geometry(communes), col="lightblue",border=NA)
typoLayer(actifs, var="typo",
legend.values.order = c("Plus d'ouvriers que de cadres",
"Plus de cadres que d'ouvriers",
"Indéterminé" ),
col = c(colouv, colsup, 'grey'),
border = "white", lwd=0.2, legend.title.txt = "Type dominant",
add = T)
labelLayer(x = actifs[data.frame(actifs)[,"total"] %in% top[1:15],], txt = "name", halo=TRUE, cex = 0.6, col= "#000000", bg = "#FFFFFF50", overlap = FALSE)
title("Ouvriers et cadres en Occitanie, 2014")A vous de jouer !
On soite réaliser une carte sur le chômage des ouvriers en Occitanie. Réalisez une carte avec les parametres suivants :
- Mode de représentation : cartograme de Dorling combinée + typologie.
- Typologie : 2 catégories (communes ou le taux de chômage des ouvriers est supérieur/inférieur à 20 %)
- Cercles proportionnels : population ouvrière (actifs + chomeurs)
- Nommer les communes où plus de 500 ouvriers sont au chômage
Voici le code à exécuter pour calculer cette typologie (même démarche que plus haut)
communes$typo <- "Indéterminé"
communes$r <- communes$ouv_chom / communes$ouv
communes$r[is.na(communes$r)] <- 0
communes[communes$r >= 0.2,"typo"] <- "Supérieur à 20 %"
communes[communes$r < 0.2,"typo"] <- "Inférieur à 20 %"Exo 3 - Cartographier l’Europe
Objectifs
- Thématique : Réaliser des cartes de densité de population et de PIB par habitant pour les régions européennes.
- Espace d’étude : Europes, NUTS3
- Visualisation associée : cartes choroplèthes et cartes en grilles
- packages utilisés : sf, cartography, eurostat, reshape2
Créez un nouveau fichier “Exo3.R”.
Chargement et préparation des données
Chargez les géométries qui seront utilisées pour visualiser vos résultats. Il s’agit des géométries de référence des territoires européens adaptés à la cartographie thématique proposés par l’UMS RIATE (ref )
# Fond de carte principal
regions <- st_read(dsn = "data/Europe/GREAT/GEOM_EU34_NUTS13/geom_eu34_nuts13.shp", stringsAsFactors = F) ## Reading layer `geom_eu34_nuts13' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/Europe/GREAT/GEOM_EU34_NUTS13/geom_eu34_nuts13.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1480 features and 7 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 2641758 ymin: 1352640 xmax: 7313157 ymax: 5411284
## epsg (SRID): NA
## proj4string: +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
# Elements d'habillage
borders <- st_read("data/Europe/GREAT/OTHERS/borders_EU34.shp", stringsAsFactors = F)## Reading layer `borders_EU34' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/Europe/GREAT/OTHERS/borders_EU34.shp' using driver `ESRI Shapefile'
## Simple feature collection with 96 features and 8 fields
## geometry type: MULTILINESTRING
## dimension: XY
## bbox: xmin: 2772939 ymin: 1743120 xmax: 5805008 ymax: 5304194
## epsg (SRID): NA
## proj4string: +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
countries <- st_read("data/Europe/GREAT/OTHERS/background_countries.shp",stringsAsFactors = F)## Reading layer `background_countries' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/Europe/GREAT/OTHERS/background_countries.shp' using driver `ESRI Shapefile'
## Simple feature collection with 73 features and 1 field
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 2600302 ymin: 1249109 xmax: 7346435 ymax: 5455936
## epsg (SRID): NA
## proj4string: +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
coastlines <- st_read("data/Europe/GREAT/OTHERS/coast_wide.shp", stringsAsFactors = F)## Reading layer `coast_wide' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/Europe/GREAT/OTHERS/coast_wide.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## geometry type: MULTILINESTRING
## dimension: XY
## bbox: xmin: 2611445 ymin: 1249569 xmax: 7345439 ymax: 5454917
## epsg (SRID): NA
## proj4string: +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
cyprus <- st_read("data/Europe/GREAT/OTHERS/cyprus_north_part.shp", stringsAsFactors = F)## Reading layer `cyprus_north_part' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/Europe/GREAT/OTHERS/cyprus_north_part.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 6380204 ymin: 1646065 xmax: 6523706 ymax: 1761964
## epsg (SRID): NA
## proj4string: +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
remote <- st_read("data/Europe/GREAT/OTHERS/remote_territories_wide.shp", stringsAsFactors = F)## Reading layer `remote_territories_wide' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/Europe/GREAT/OTHERS/remote_territories_wide.shp' using driver `ESRI Shapefile'
## Simple feature collection with 2 features and 1 field
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 2657347 ymin: 3613414 xmax: 2881832 ymax: 3836004
## epsg (SRID): NA
## proj4string: +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
seaboxes <- st_read ("data/Europe/GREAT/OTHERS/remote_area_background_wide.shp",stringsAsFactors = F)## Reading layer `remote_area_background_wide' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/Europe/GREAT/OTHERS/remote_area_background_wide.shp' using driver `ESRI Shapefile'
## Simple feature collection with 8 features and 23 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 2657347 ymin: 2816552 xmax: 2881832 ymax: 4632866
## epsg (SRID): NA
## proj4string: +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
background <- st_read("data/Europe/GREAT/OTHERS/background.shp", stringsAsFactors = F)## Reading layer `background' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/Europe/GREAT/OTHERS/background.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 2600302 ymin: 1249109 xmax: 7346435 ymax: 5455936
## epsg (SRID): NA
## proj4string: +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
Il est toujours utile de visualiser les géométries utilisées et consulter la table attributaire avant de lancer les traitements…
# Ajuster les marges
par(mar = c(0.5,0.5,1.5,0.5))
# Visualisation des couches géographiques
plot(st_geometry(background), col = "#e0ecff", border = NA, add =FALSE)
plot(st_geometry(countries), col = "#e6e6e6", border = NA, add =TRUE)
plot(st_geometry(seaboxes), col = "#f7fcfe", border = NA, add = TRUE)
plot(st_geometry(regions), col = "#7ccdea", border = "#f7fcfe", lwd = 0.15, add=T) # Les representations cartographiques porteront sur cette couche
plot (st_geometry(remote), col = "#e6e6e6", border =NA, add=TRUE)
plot (st_geometry(borders), col = "#f7fcfe",lwd = 0.5, add = TRUE)
plot (st_geometry(cyprus), col = "#f7fcfe", border = NA, add = TRUE)
plot (st_geometry(coastlines), col = "#d2dbef", lwd = 0.15, add = TRUE)
plot(st_geometry(seaboxes), col = NA, border = "#bbbdc0", lwd = 0.15, add = TRUE)# Table attributaire de la couche cible pour les traitements cartographiques
head(regions)## Simple feature collection with 6 features and 7 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 4655611 ymin: 2655076 xmax: 4855527 ymax: 2823072
## epsg (SRID): NA
## proj4string: +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
## id name level0 level1 level2 level3 level23
## 1 AT111 Mittelburgenland 0 0 0 1 0
## 2 AT112 Nordburgenland 0 0 0 1 0
## 3 AT113 Südburgenland 0 0 0 1 0
## 4 AT121 Mostviertel-Eisenwurzen 0 0 0 1 0
## 5 AT122 Niederösterreich-Süd 0 0 0 1 0
## 6 AT123 Sankt Pölten 0 0 0 1 0
## geometry
## 1 MULTIPOLYGON (((4819333 274...
## 2 MULTIPOLYGON (((4819333 274...
## 3 MULTIPOLYGON (((4808360 271...
## 4 MULTIPOLYGON (((4729815 281...
## 5 MULTIPOLYGON (((4762433 279...
## 6 MULTIPOLYGON (((4768449 279...
Ce bloc de code est dédié est à la recolte des données disponibles sur le site d’Eurostat grâce au package eurostat. La procédure à suivre pour télécharger ses propres données est la suivante :
- Consultez le contenu de la base de données Eurostat et identifiez la table qui contient vos indicateurs d’intérêt (cf ci-dessous)
- Importez la librairie eurostat.
- Téléchargez la table d’intérêt avec la fonction get_eurostat. C’est un tibble.
- Consultez les dimensions du tibble avec la fonction str()
- Filtrez les dimensions du tibble afin de réduire sa structure à deux dimensions (objets géographiques * indicateur par exemple).
- Transformez le tibble en dataframe avec la fonction dcast(). Dans notre cas de figure nous conservons en ligne les objets géographiques et en colonne la dimension temporelle de l’indicateur d’intérêt.
- Quand tous les indicateurs ont été téléchargés et préparés, les joindre à vos géométries de référence (regions dans notre cas de figure).
library("eurostat")
library("reshape2")
# POPULATION
var <- "demo_r_pjangrp3" # Table Eurostat d'intérêt
data <- get_eurostat(var, time_format = "num") # Telecharger la table ESTAT
data <- subset(data, data$sex == "T") # Filtre des dimensions du tibble
data <- subset(data, data$age == "TOTAL") # Filtre des dimensions du tibble
data <- dcast(data, geo ~ time, value.var = "values") # Transformation en dataframe des dimensions restantes
colnames(data) <- c("geo","POP_2014","POP_2015","POP_2016","POP_2017") # Renommer les colonnes du dataframe de façon explicite
# SURFACE
var <- "reg_area3"
area <- get_eurostat(var, time_format = "num") # Telecharger la table ESTAT
area <- subset(area, area$landuse == "TOTAL") # Filtre des dimensions de la table
area <- dcast(area, geo ~ time, value.var = "values")
colnames(area) <- c("geo","AREA")
# PIB
var <- "nama_10r_3gdp"
gdp <- get_eurostat(var, time_format = "num")
gdp <- subset(gdp, gdp$unit == "MIO_EUR")
gdp <- dcast(gdp, geo ~ time, value.var = "values")
fields <- c("geo", "2014", "2015", "2016") # On ne garde que les valeurs de PIB correspondant aux valeurs de population disoponibles
gdp <- gdp[,fields]
colnames(gdp) <- c("geo","GDP_2014","GDP_2015","GDP_2016")
# Jointure
regions <- merge(x = regions, y = data, by.x = "id", by.y = "geo", all.x=T) # population
regions <- merge(x = regions, y = area, by.x = "id", by.y = "geo", all.x=T) # surface
regions <- merge(x = regions, y = gdp, by.x = "id", by.y = "geo", all.x=T) # surfaceImportez le jeu de données créé suite à cette procédure, associé aux géométries européennes de réference.
regions <- st_read("data/Europe/regions_data.shp", stringsAsFactors = F)## Reading layer `regions_data' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/Europe/regions_data.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1480 features and 15 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 2641758 ymin: 1352640 xmax: 7313157 ymax: 5411284
## epsg (SRID): NA
## proj4string: +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
Les données sont prêtes pour la visualisation !
Pour choisir l’année de référence sur laquelle portera la visualisation, regardez le nombre de valeurs manquantes…
apply(regions,2,function(x) sum(is.na(x)))## id name level0 level1 level2 level3 level23 POP_2014
## 0 0 0 0 0 0 0 2
## POP_2015 POP_2016 POP_2017 AREA GDP_2014 GDP_2015 GDP_2016 geometry
## 1 1 1 1 111 113 1161 0
… Et calculez le ratio sur lequel va porter la visualisation en conséquence
regions$DENS_2017 <- regions$POP_2017 / regions$AREAVisualisation : cartes choroplèthes
La variable visuelle adaptée pour représenter des données quantitatives relatives est la couleur ordonnée (ou niveaux de gris). Le traitement de ce type de données implique au prélable de définir le nombre de classes, une méthode de discrétisation et une palette de couleurs sur laquelle portera la représentation.
La discrétisation consiste à découper une série statistique en classes de valeurs. Il n’existe pas de méthode optimum. Chaque méthode donnera une carte différente, plus ou moins conforme à la distribution de départ. L’agrégation de données en classes introduit de fait une erreur ou une distorsion dans la perception de cette distribution. Le choix d’une méthode de discrétisation dépend donc des propriétés de la distribution mais aussi des objectifs cartographiques que l’on s’est fixés.
Dans cartography, plusieurs méthodes de discrétisation sont proposées. Jetez un oeil à la documentation de la fonction getBreaks pour en savoir plus.
On peut maintenant observer avec un graphique simple les effets de plusieurs méthodes de discrétisation sur notre variable de densité de population.
Compte-tenu du nombre d’unités territoriales NUTS3 en Europe, nous choisissons 10 classes pour guider la discrétisation.
library(cartography)
# 4 graphiques par plot
par(mfrow=c(2,2))
# Methode discretisation (quantiles)
breaks <- getBreaks(v = regions$DENS_2017, nclass = 10, method = "quantile")
hist(regions$DENS_2017, probability = TRUE, breaks = breaks, col = "#F0D9F9", xlim = c(0,4000), main = "Quantiles")
abline(v = median(regions$DENS_2017, na.rm = T), col = "blue", lwd = 3)
abline(v = mean(regions$DENS_2017, na.rm = T), col = "red", lwd = 3)
# Methode discretisation (intervalles égaux)
breaks <- getBreaks(v = regions$DENS_2017, nclass = 10, method = "equal")
hist(regions$DENS_2017, probability = TRUE, breaks = breaks, col = "#F0D9F9", xlim = c(0,4000), main = "Intervalles égaux")
abline(v = median(regions$DENS_2017, na.rm = T), col = "blue", lwd = 3)
abline(v = mean(regions$DENS_2017, na.rm = T), col = "red", lwd = 3)
# Methode discretisation (moyenne-écart type)
breaks <- getBreaks(v = regions$DENS_2017, nclass = 10, method = "msd", k=1, middle = TRUE)
hist(regions$DENS_2017, probability = TRUE, breaks = breaks, col = "#F0D9F9", xlim = c(0,4000), main = "Moyenne/écart-type")
med <- median(regions$DENS_2017, na.rm = T)
abline(v = median(regions$DENS_2017, na.rm = T), col = "blue", lwd = 3)
abline(v = mean(regions$DENS_2017, na.rm = T), col = "red", lwd = 3)
# Methode discretisation (geom)
breaks <- getBreaks(v = regions$DENS_2017, nclass = 10, method = "geom")
hist(regions$DENS_2017, probability = TRUE, breaks = breaks, col = "#F0D9F9", xlim = c(0,4000), main = "Progression géométrique")
abline(v = median(regions$DENS_2017, na.rm = T), col = "blue", lwd = 3)
abline(v = mean(regions$DENS_2017, na.rm = T), col = "red", lwd = 3)Explorez la documentation de la fonction carto.pal de cartography pour choisir votre palette. Utilisez l’instruction display.carto.all pour visualiser les palettes disponibles. Nous choisissons ici la palette “red.pal”
# Choix de la palette
display.carto.all(10)cols <- carto.pal(pal1 = "red.pal", n1 = 10)Tout est prêt ! Cartographions cet indicateur avec la méthode des quantiles, la palette “red” en 10 classes. Associons des couches d’habillage. Et mettons en page la carte (flèche du nord, échelle, sources) grâce à l’instruction layoutLayer de cartography.
# Ajuster les marges et les paramètes de visu
par(mar = c(0.5,0.5,1.5,0.5))
par(mfrow=c(1,1))
# Mise en page (back)
plot(st_geometry(background), col = "#e0ecff", border = NA, add =FALSE)
plot(st_geometry(countries), col = "#c6c4c4", border = NA, add =TRUE)
plot(st_geometry(seaboxes), col = "#f7fcfe", border = NA, add = TRUE)
# Cartographie (ratio)
choroLayer(x = regions,
var = "DENS_2017",
method ="quantile",
nclass = 10,
col = cols,
border = NA,
add = TRUE,
legend.pos = "topleft",
legend.title.txt = "Densité de population, 2017 (hab/km²)",
legend.values.rnd = 0)
# Mise en page (top)
plot (st_geometry(remote), col = "#e6e6e6", border =NA, add=TRUE)
plot (st_geometry(borders), col = "#f7fcfe",lwd = 0.5, add = TRUE)
plot (st_geometry(cyprus), col = "#f7fcfe", border = NA, add = TRUE)
plot (st_geometry(coastlines), col = "#d2dbef", lwd = 0.15, add = TRUE)
plot(st_geometry(seaboxes), col = NA, border = "#bbbdc0", lwd = 0.15, add = TRUE)
# Sources + elements d'habillage
layoutLayer(title = "Densité de population dans les NUTS3 européens, 2017",
author = "Réalisation : Serial Mappers, 2018", sources = "Source, Eurostat, 2018",
scale = 200, tabtitle = TRUE,
frame = FALSE,theme = "red.pal",
south = TRUE)A vous de jouer !
Reproduisez cette carte avec les paramètres suivants :
- Discrétisation : progression géométrique
- Nombre de classes : 8
- Palette de couleurs : Bleue
Visualisation : cartes carroyées
La méthode du carroyage consiste à découper l’espace géographique en un maillage formé de carrés réguliers dans une projection donnée. Elle permet de s’affranchir de l’arbitraire et de l’irrégularité d’un découpage administratif, comme c’est ici le cas avec l’Europe au niveau NUTS3. La donnée est répartie sur ce quadrillage régulier au prorate de la surface représentée.
Comment faire ? L’instruction getGridLayer du package cartography permet de construire ces grilles régulières. Expolorez la documentation pour vous sensibiliser aux arguments de la fonction.
La procédure est la suivante :
- Définir le pas de grille (en mètres), son type et les variables qui sont transformées
- Pour les calculs de densité, convertir la surface en km².
- Pour les autres indicateurs, combiner les données de la façon suivante : ratio = stock1/stock2
- Choisir sa palette, son nombre de classes et sa méthode de discrétisation
- Représenter la grille avec la fonction choroLayer de cartography.
# Transformation des donnees dans une grille reguliere de 50 km
mygrid50k <- getGridLayer(x = regions, cellsize = 50000 * 50000,
type = "regular", var = c("POP_2017","POP_2014","GDP_2014"))
# Transformation des donnees dans une grille reguliere de 100 km
mygrid100k <- getGridLayer(x = regions, cellsize = 100000 * 100000,
type = "regular", var = c("POP_2017","POP_2014","GDP_2014"))
# Conversion en km²
mygrid50k$densitykm <- mygrid50k$POP_2017 * 1000 * 1000 / mygrid50k$gridarea
mygrid100k$densitykm <- mygrid100k$POP_2017 * 1000 * 1000 / mygrid100k$gridarea
# Cartographie
par(mfrow=c(1,2)) # 2 cartes
par(mar = c(0.5,0.5,1.5,0.5)) # Ahuster les marges
cols <- carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 5, n2 = 5)
choroLayer(x = mygrid50k, var = "densitykm",
border = "grey80",col = cols, method ="quantile",
nclass = 10, legend.pos = "topright", legend.values.rnd = 1,
legend.title.txt = "Densité de population, 2017")
layoutLayer(title = "Densité de population dans une grille régulière de 50km",
author = "Serial Mappers, 2018", sources = "Eurostat, 2018",
scale = 200,
frame = TRUE,
col = "red",
coltitle = "white",
south = TRUE)
choroLayer(x = mygrid100k, var = "densitykm",
border = "grey80",col = cols, method ="quantile",
nclass = 10, legend.pos = "topright", legend.values.rnd = 1,
legend.title.txt = "Densité de population, 2017")
layoutLayer(title = "Densité de population dans une grille régulière de 100km",
author = "Serial Mappers, 2018", sources = "Eurostat, 2018",
scale = 200,
frame = TRUE,
col = "red",
coltitle = "white",
south = TRUE)A vous de jouer !
Créez une nouvelle carte carroyée avec les paramètres suivants :
- Indicateur : PIB / habitant 2014
- Résolution de grille : 100 km
- Type de grille : hexagones (nb : il faut charger la librairie sp au préalable)
- Discrétisation : quantiles (8 classes)
- Palette de couleurs : Double palette (vert/rouge)
Exécutez ce code qui permet au prélable de préparer la grille et retirer les valeurs manquantes qui nuiraient au calcul.
library(sp)
# Retirer les valeurs manquantes
mygrid100k <- regions
mygrid100k <- mygrid100k[!is.na(mygrid100k$GDP_2014),]
mygrid100k <- mygrid100k[!is.na(mygrid100k$POP_2014),]
# Creation de la grille
mygrid100k <- getGridLayer(x = mygrid100k, cellsize = 100000 * 100000,
type = "hexagonal", var = c("POP_2014","GDP_2014"))
# Calcul PIB par habitant
mygrid100k$GDP_HAB <- mygrid100k$GDP_2014 * 1000000 / mygrid100k$POP_2014Exo 4 - Cartographier le Monde
Objectifs :
- Thématique : Les inégalités d’IDH et/ou d’espérance de vie
- Espace d’étude : Le Monde, pays.
- Visualisation associée : carte chorplèthe, discontinuités
- packages utilisés : sf, cartography, rnaturalearth, readxl, countrycode
Import du fond de carte
méthode 1 : à la main (après téléchargement sur le site Natural Earth)
library("sf")
countries <- st_read(dsn = "data/world/naturalearth/ne_110m_admin_0_countries.shp", stringsAsFactors = F) Méthode 2 : directement via R
url <- "https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/cultural/ne_110m_admin_0_countries.zip"
download.file(url = url, destfile = "data/download/countries.zip")
unzip("data/download/countries.zip", exdir = "data/download", unzip = "internal")
file.remove("data/download/countries.zip")
countries <- st_read(dsn = "data/download/ne_110m_admin_0_countries.shp", stringsAsFactors = F)Méthode 3 : directement via R
library("rnaturalearth")
countries <- ne_countries(scale = 50, type = "countries", continent = NULL,
country = NULL, geounit = NULL, sovereignty = NULL,
returnclass = "sf")Le fond de carte peut être visualisé avec l’instruction plot
plot(st_geometry(countries))On ne conserve que les champs utiles. On les renomme.
countries <- countries[,c("adm0_a3", "admin", "geometry")]
colnames(countries) <- c("id","name","geometry")On Supprime l’Antarctique
countries <- countries[countries$id != "ATA",]On supprime les lignes vides
countries <- countries[!is.na(countries$id),]plot(st_geometry(countries))Import des données atributaires
Import des données sur l’éspérance de vie (source : Banque mondiale)
library("readxl")
file <-"data/world/worldbank/API_SP.DYN.LE00.IN_DS2_en_excel_v2_10081535.xls"
sheet <- "Data"
lifexp <- data.frame(read_excel(file, skip = 3, sheet = sheet))
lifexp <- lifexp[,c("Country.Code","Country.Name","X2016")]
colnames(lifexp) <- c("id","name","lifexp2016")
lifexp$lifexp2016 <- as.numeric(as.character(lifexp$lifexp2016))Import des données sur l’Indice de deveoppement humain (source : Human Development Report)
file <-"data/world/hdr/Human Development Index (HDI).csv"
hdi <- read.csv2(file = file, sep = ",",skip = 1)
hdi <- hdi[,c("Country","X2016")]head(hdi,4)## Country X2016
## 1 Afghanistan 0.494
## 2 Albania 0.782
## 3 Algeria 0.753
## 4 Andorra 0.856
Le tableau de données sur l’IDH ne possède pas de codes géographiques, mais uniquement le nom des pays. Cela rend la jointure avec le fond de carte compliquée. Une solution : le package countrycode
library("countrycode")
hdi$id <- countrycode(hdi$Country, 'country.name', 'iso3c')## Warning in countrycode(hdi$Country, "country.name", "iso3c"): Some values were not matched unambiguously: Eswatini (Kingdom of)
Un peu de mise en forme. Les codes sont bien là. C’est presque magique… Mais attention, dans la pratique, une vérification manuelle s’impose.
hdi <- hdi[,c("id","Country","X2016")]
colnames(hdi) <- c("id","name","hdi2016")
hdi$hdi2016 <- as.numeric(as.character(hdi$hdi2016))
head(hdi,4)## id name hdi2016
## 1 AFG Afghanistan 0.494
## 2 ALB Albania 0.782
## 3 DZA Algeria 0.753
## 4 AND Andorra 0.856
Cartographie
Chargement du package
library("cartography")Changer la projection du fond de carte (voir le site spatialreference
countries <- st_transform(countries, "+proj=cea +lon_0=0 +lat_ts=45 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs ")
plot(st_geometry(countries))Jointures des deux fichiers de données
countries <- merge(countries, lifexp, by="id", all.x=TRUE)
countries <- merge(countries, hdi, by="id", all.x=TRUE)
countries <- countries[,c("id","name.x","lifexp2016","hdi2016","geometry")]
colnames(countries)[2] <- "name"On travaille sur l’indice de développement humain
analyse de la distribution
var <- countries$hdi2016
hist(countries$hdi2016, probability = TRUE, nclass = 30)On choisit une discretisation par la méthode des quantiles sans classe centrale
breaks <- getBreaks(v = var, nclass = 6, method = "quantile")
cols <- carto.pal(pal1="blue.pal", n1 = 3, pal2 = "orange.pal", n2 = 3, middle = FALSE, transparency = FALSE)Réalisation de la carte
layoutLayer(extent = countries,
bg="#EEEEEE",
title = "Indice de développement humain en 2016",
scale = NULL,
sources = "Human development Report, 2018",
author = "les serial mappers"
)
choroLayer(x = countries,
var = "hdi2016",
breaks = breaks,
col = cols,
border = "white",
lwd=0.2,
add = TRUE,
legend.pos = "topleft",
legend.title.txt = "IDH, 2016",
legend.values.rnd = 2)Et si on ajoutait des lignes de discontinuités à cette carte?
L’instruction getBorders (du package cartography) permet d’extraire les frontières entre les unités spatiales.
borders <- getBorders(countries)
plot(st_geometry(borders), col="black", lwd=1)Quid des discontinuités maritimes ? Par exemple entre l’Italie et la Tunisie ? Ou entre le Maroc et l’Espagne ? L’instruction getOuterBorders (package cartography) permet de les générer. Cette instruction peut prendre un peu de temps.
outer <- getOuterBorders(x = countries, res = 100000, width=500000)
plot(st_geometry(borders), col="#CCCCCC", lwd=1)
plot(st_geometry(outer), col="red", lwd=1, add=T)On assemble les deux couches avec rbind
b <- rbind(borders,outer)Tout est prêt pour réaliser une carte de discontinuités
# 1 - la carte choroplèthe
layoutLayer(extent = countries,
bg="#EEEEEE",
title = "Indice de développement humain en 2016",
scale = NULL,
sources = "Human development Report, 2018",
author = "les serial mappers"
)
choroLayer(x = countries,
var = "hdi2016",
breaks = breaks,
col = cols,
border = "white",
lwd=0.2,
add = TRUE,
legend.pos = "topleft",
legend.title.txt = "IDH, 2016",
legend.values.rnd = 2)
# 2 - les discontinuités
discLayer(x = b, df = countries,
var = "hdi2016", col="black", nclass=4,
method="quantile", threshold = 0.25, sizemin = 0.2,
sizemax = 8, type = "abs",
legend.title.txt = "Discontinuities\nabsolues",
legend.pos = "bottomleft", add=TRUE)BONUS
Utiliser des données OSM
library("sf")
library("cartography")
depts <- st_read(dsn = "data/France/IGN/DEPARTEMENT.shp", stringsAsFactors = F)## Reading layer `DEPARTEMENT' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/France/IGN/DEPARTEMENT.shp' using driver `ESRI Shapefile'
## Simple feature collection with 96 features and 11 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 99217.1 ymin: 6049646 xmax: 1242417 ymax: 7110480
## epsg (SRID): NA
## proj4string: +proj=lcc +lat_1=44 +lat_2=49 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
colnames(depts)## [1] "ID_GEOFLA" "CODE_DEPT" "NOM_DEPT" "CODE_CHF" "NOM_CHF"
## [6] "X_CHF_LIEU" "Y_CHF_LIEU" "X_CENTROID" "Y_CENTROID" "CODE_REG"
## [11] "NOM_REG" "geometry"
Paris <- depts[depts$CODE_DEPT == 75,]
plot(st_geometry(Paris))Paname <- getTiles(x = Paris, type = "osm", crop = TRUE)
Paname <- getTiles(x = Paris, type ="opencycle", crop = TRUE)
tilesLayer(Paname)
plot(st_geometry(Paris), lwd=4, add=T)library("osmdata")## Data (c) OpenStreetMap contributors, ODbL 1.0. http://www.openstreetmap.org/copyright
prj <- st_crs(Paris)
bbox <- st_bbox(st_transform(Paris,4326))
q <- opq(bbox = bbox , timeout = 2000) %>% add_osm_feature(key = 'man_made', value = 'surveillance')
cameras <- osmdata_sf(q)$osm_points
cameras <- st_transform(cameras, prj)
cameras$ok <- st_intersects(st_geometry(cameras), st_geometry(Paris), sparse = FALSE)
cameras <- cameras[cameras$ok == T,]
Paname <- getTiles(x = Paris, type ="cartodark", crop = TRUE)
tilesLayer(Paname)
plot(st_geometry(Paris), lwd=1, border="white", lty=2, add=T)
plot(st_geometry(cameras), pch=20, col="red", add=T)Généraliser un fond de carte
Bla bla sur la généralisation cartographique
library("rmapshaper")Import et affichage d’un fond de carte
countries <- st_read(dsn = "data/world/naturalearth/ne_110m_admin_0_countries.shp", stringsAsFactors = F)## Reading layer `ne_110m_admin_0_countries' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/world/naturalearth/ne_110m_admin_0_countries.shp' using driver `ESRI Shapefile'
## Simple feature collection with 177 features and 94 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
plot(st_geometry(countries), col="green")Généralisation avec sf (ca marche mal)
countries2 <- st_simplify(x = countries, preserveTopology = TRUE, dTolerance = 2)## Warning in st_simplify.sfc(st_geometry(x), preserveTopology, dTolerance):
## st_simplify does not correctly simplify longitude/latitude data, dTolerance
## needs to be in decimal degrees
plot(st_geometry(countries2), col="green")Généralisation avec rmapshaper (Youpie !!!)
countries2 <- ms_simplify(countries, keep = 0.1)
plot(st_geometry(countries2), col="green")BIBLIO
- BEGUIN Michèle, PUMAIN Denise, La représentation des données géographiques, Statistique et cartographie, coll. Cursus, Armand Colin, nouvelle édition 2000, 192p.
- BERTIN Jacques, Sémiologie graphique, Monton-Gauthier-Villars, 1967, 1973, 432 p. Disponible à ce jour dans la collection Réimpression de l’EHESS, 1998.
- BERTIN Jacques, La graphique et le traitement graphique de l’information, Flammarion, 1977, 280 p.
- BRUNET Roger, La carte mode d’emploi, Fayard-Reclus, 1987, 270p.
- ElementR (Groupe), R et espace, traitement de l’information géographique, 2014, Framabook.
- LAMBERT Nicolas : Carnet de recherche neocarto
- LAMBERT Nicolas, ZANIN Christine : « Manuel de cartographie. Principes, méthodes, applications », cursus, Armand Colin, 18 mai 2016, 224p.
- Giraud Timothée : Carnet de recherche rgeomatic
- ZANIN Christine, TREMELO Marie-Laure, Savoir faire une carte, Aide à la conception et à la réalisation d’une carte thématique univariée, coll. Sup géographie, Belin, 2003, 200p.
— R infos —
Version de R
R.version## _
## platform x86_64-pc-linux-gnu
## arch x86_64
## os linux-gnu
## system x86_64, linux-gnu
## status
## major 3
## minor 4.4
## year 2018
## month 03
## day 15
## svn rev 74408
## language R
## version.string R version 3.4.4 (2018-03-15)
## nickname Someone to Lean On
Session
sessionInfo()## R version 3.4.4 (2018-03-15)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.1 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=fr_FR.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=fr_FR.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=fr_FR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] rmapshaper_0.4.0 osmdata_0.0.7.001 countrycode_1.00.0
## [4] rnaturalearth_0.1.0 sp_1.3-1 cartogram_0.1.0
## [7] cartography_2.1.2 readxl_1.1.0 sf_0.6-3
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.18 lubridate_1.7.4
## [3] lattice_0.20-35 png_0.1-7
## [5] class_7.3-14 rprojroot_1.3-2
## [7] digest_0.6.17 V8_1.5
## [9] mime_0.5 R6_2.2.2
## [11] cellranger_1.1.0 plyr_1.8.4
## [13] backports_1.1.2 rnaturalearthdata_0.1.0
## [15] evaluate_0.11 e1071_1.7-0
## [17] geojson_0.2.0 httr_1.3.1
## [19] highr_0.7 pillar_1.3.0
## [21] rlang_0.2.2 lazyeval_0.2.1
## [23] curl_3.2 geojsonlint_0.2.0
## [25] rstudioapi_0.7 miniUI_0.1.1.1
## [27] raster_2.6-7 geojsonio_0.6.0
## [29] rmarkdown_1.10 rgdal_1.3-4
## [31] foreign_0.8-69 jqr_1.0.0
## [33] stringr_1.3.1 questionr_0.6.3
## [35] shiny_1.1.0 rosm_0.2.2
## [37] compiler_3.4.4 httpuv_1.4.5
## [39] xfun_0.3 rgeos_0.3-28
## [41] htmltools_0.3.6 tibble_1.4.2
## [43] bookdown_0.7 jsonvalidate_1.0.0
## [45] crayon_1.3.4 later_0.7.5
## [47] grid_3.4.4 spData_0.2.9.4
## [49] jsonlite_1.5 xtable_1.8-3
## [51] DBI_1.0.0 magrittr_1.5
## [53] units_0.6-1 rmdformats_0.3.3
## [55] stringi_1.2.4 promises_1.0.1
## [57] xml2_1.2.0 packcircles_0.3.3
## [59] tools_3.4.4 abind_1.4-5
## [61] yaml_2.2.0 maptools_0.9-4
## [63] classInt_0.2-3 rvest_0.3.2
## [65] knitr_1.20